Infinite disorder scaling of random quantum magnets in three and higher dimensions 
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Using a very efficient numerical algorithm of the strong disorder renormalization group method 
we have extended the investigations about the critical behavior of the random transverse-field Ising 
model in three and four dimensions, as well as for Erdos-Renyi random graphs, which represent 
infinite dimensional lattices. In all studied cases an infinite disorder quantum critical point is 
identified, which ensures that the applied method is asymptotically correct and the calculated critical 
exponents tend to the exact values for large scales. We have found that the critical exponents are 
independent of the form of (ferromagnetic) disorder and they vary smoothly with the dimensionality. 



I. INTRODUCTION 

Quantum phase transitions are among the fundamen- 
tal problems of modern physics, the properties of which 
are studied in solid state physics, quantum field-theory, 
quantum information and statistical mechanics^. These 
transitions take place at T = temperature, i.e. in 
the ground state of the quantum system by varying a 
control parameter, such as the strength of a transverse 
field. One basic question in this field of research is 
how quenched disorder influences the properties of quan- 
tum phases and phase transitions. In this respect quan- 
tum spin glasses and the glass transition are particularly 
interesting 2 -. This latter problem theoretically is very 
challenging, since the corresponding quantum state is the 
result of an interplay between quantum and disorder fluc- 
tuations, strong correlations and frustration. 

One of the paradigmatic models of random quan- 
tum magnets with a discrete symmetry is the random 
transverse-field Ising model (RTIM) , which is defined by 
the Hamiltonian: 
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Here the a*' z are Pauli-matrices and i, j denote sites of a 
lattice (or a graph). Experimentally the RTIM is closely 
related to the compound^ LiHo x Yi_ x F4, in which there 
is a dipolc-coupling between the Ising spins, thus the in- 
teraction is long-ranged. Applying a magnetic field H t 
transverse to the Ising axis results in a transverse field 
of strength, hi = Hf, but this transverse field induces a 
random longitudinal field^ via the off-diagonal terms of 
the dipolar interaction. In the theoretical investigations 
the interactions in the RTIM arc generally assumed to be 
short-ranged, thus the first sum in Eq.([T]) runs over near- 
est neighbors. Furthermore the couplings and the hi 
transverse fields are independent random numbers, which 
are taken from the distributions, p(J) and q(h), respec- 
tively. For random ferromagnets we have J > 0, whereas 
for spin-glasses there are both ferro- and antiferromag- 
netic couplings. Here we are basically interested in the 



former problem. 

Detailed theoretical results about the RTIM are known 
in one dimension (ID) due to a complete analytical so- 
lution of a renormalization group (RG) treatment. The 
RG results are expected to be asymptotically exact in 
the vicinity of the critical point (and also in the Griffiths- 
phase, as long as dynamical singularities are concerned^), 
which is indeed demonstrated by a comparison with in- 
dependent analytical^ and numerical*^ works. One 
important observation, that the critical properties of the 
ID model are governed by an infinite disorder fixed point 
(IDFP), in which the strength of disorder growths with- 
out limit during renormalizationii and thus become dom- 
inant over quantum fluctuations. 

The IDFP scenario is found to be valid for the 2D 
RTIM, too, as observed in numerical RG studies^— 
and in Monte Carlo (MC) simulations^. The calculated 
critical exponents are in agreement with the MC results 
about the 2D random contact process 2 ^, which is a sim- 
ple noncquilibrium model of spreading infections. The 
d-dimensional random contact process is expected to be 
in the same universality class2i as the RTIM, at least for 
strong enough disorder. 

In three dimensions, which is connected to real quan- 
tum magnets, no quantitative results are known, so far. 
Analysis of the numerical RG trajectories lead to the 
conclusion^ 2 -, that the critical behavior in this case is 
probably controlled by an IDFP, but no estimates about 
the critical exponents are available. For even higher di- 
mensions it is unclear, if the IDFP scenario stays valid 
for any finite value of d, or there is some upper critical 
dimension, d c , so that for d> d c the critical behavior is 
of conventional disorder type. We note that the largc-G? 
limit of the problem is qualitatively relevant for models 
with long-range interactions. In this respect the critical 
behavior of the (non-random) L1H0F4 system has been 
the subject of intensive MC simulations 2 ^. 

In this paper we extend the investigations about the 
critical behavior of the RTIM into the hitherto unex- 
plored three and higher dimensions. Here we have devel- 
oped a considerably improved numerical algorithm of the 
strong disorder RG (SDRG) procedure and study large 
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samples with N > 10 6 sites. In 3D and 4D we consider 
hypercubic lattices and we study the large-D limit of the 
problem, too. This latter is realized by Erdos-Rcnyi (ER) 
random graphs^ 3 - consisting of N sites and kN/2 edges 
(k > 2), which are at random positions. The sizes of 
the largest systems we studied are shown in Table HI We 
note that the numerical algorithm of the SDRG method 
is completely different from that used in the 2D case, 
however, the steps used to calculate the critical param- 
eters as well as the method of analysing the results are 
similar to that used in 2Di&. 

The structure of our paper is the following. The SDRG 
method and its improved algorithm used in this paper is 
described in SecHU Results about the critical parameters 
are calculated in Sec lIIII and discussed in Sec lIVI 



II. SDRG PROCEDURE 

In the calculation we used the SDRG procedure^, 
which has been introduced by Ma, Dasgupta and Hv±2£. 
In this method, which works for random fcrromagnets, at 
each step of the renormalization the largest local term in 
the Hamiltonan (either a coupling or a transverse field) 
is eliminated and new terms are generated between re- 
maining sites by second-order perturbation method. Af- 
ter decimating a strong coupling, say Jy, the two con- 
nected spins form a spin cluster having an additive mo- 
ment, jl = ^ + Hj, which is placed in an effective trans- 
verse field of strength: h = 1 J . After decimating a 

Jij 

large transverse field, hi, the actual spin is eliminated 
and new effective couplings are generated between each 
pair of spins being nearest neighbors to the decimated 

site, say j and k, having a value: Jjk = Jl . tk ■ If at one 

hi 

step two parallel couplings appear between two neighbor- 
ing sites the maximum of them is taken. Application of 
this "maximum rule" is exact at an IDFP and results in 
simplifications of the RG procedure. 

Here we have developed an optimized algorithm, which 
needs t ~ 0(N log N + E) time to renormalize a cluster 
with A* sites and E edges up to the last spin, irrespective 
of the dimension and topology of the cluster—. In this 
algorithm terms in the Hamiltonian are decimated in de- 
scending order in energy and we have applied the follow- 
ing theorem for transverse field decimation. According 
to this theorem for a decimated site, i, there is always 
one relevant neighboring site, j, so that after decimating 
i only those renormalizcd couplings should be created, 
which connect j with its new neighboring sites2£. All the 
couplings which start from other nearest neighbors of i 
(and does not end at j) are irrelevant and need not be 
created. In this way during one RG step not only the 
number of sites is reduced (by one), but the number of 
couplings, as well. Using this algorithm we avoid to gen- 
erate almost fully connected clusters, which is the main 
drawback of the naive implementation of the method in 



higher dimensions^, having a performance: t ~ 0(N 3 ). 
Renormalization with the naive and the improved algo- 
rithms is illustrated in FigjTJ 




FIG. 1: (Color online) Illustration of the renormalization 
procedure. In the original cluster (left panel) the central site, 
i, is to be decimated. In the naive implementation (middle 
panel) several new couplings are generated. In the improved 
algorithm (right panel) just the couplings to the relevant site, 
j, are created. 



III. CALCULATION OF CRITICAL 
PARAMETERS 

In the actual calculation - in order to check univer- 
sality and to control the disorder dependence of the es- 
timates - we have used two different forms of random- 
ness. Both have the same uniform distribution of the 
couplings: p(J) = 8(J)0(1 — J) (Q(x) being the Heav- 
iside step- function), which are ferromagnetic. For the 
'box-/i' disorder also the transverse fields are uniformly 

distributed: q(h) = — Q(h)Q(hb — h), whereas for the 
hb 

'fixed-/i' model we have a constant transverse field^: 
hi = hf. We used the logarithmic variable, 8 = \n(hb) or 
9 = ln(/i/), as a quantum control parameter. We have 
checked that the computational time to renormalize an 
N = 10 3 (N = 10 6 ) cluster is typically - 0.015 (~ 50) 
second (in a 2.4GHz processor), which does not depend 
on the dimension and the topology of the cluster. The 
numbers of realizations used in the calculations were typ- 
ically 40000 but even for the largest sizes we have at least 
10000 samples. 

In the first step of the calculation for each random 
sample, a, we have determined a pseudo-critical point, 
6 c (a,N), by a variant of the doubling method. In this 
procedure^ we glue together two identical copies (a, a') 
of the sample by surface couplings^ and renormalize it 
up to the last site for different values of the control pa- 
rameter, 6. The renormalization is found to be qualita- 
tively different for 8 < c (a,N) and for 8 > 6 c (a,N). 
For weak quantum fluctuations, 8 < 8 c (a,N), the last 
decimated spin cluster contains equivalent sites of a and 
a' . These sites and thus the two replicas are correlated 
and we call this cluster as a correlation cluster. The mo- 
ment of the correlation cluster, fj,(a,N), goes to zero as 
the pseudo-critical point is approached. On the contrary 
for 9 > 8 c (a, N) in the last decimated spin cluster there 
are no equivalent sites of a and a' and thus there is no 
correlation cluster. 



3 



10 u 
10" 1 
10" 2 
10" 3 

10- 4 








f 








j 


t 




X 

% 




\ 



-0.4 -0.2 



0.2 0.4 
L=16 + 
24 

32 * 
48 □ 
64 

96 o 



-0.1 



-0.09 -0.08 



-0.07 



-0.06 -0.05 -0.04 



FIG. 2: (Color online) Distribution of the pseudo-critical 
points, 8 c {a,N), for various sizes, N — L 3 , for fixed-ft ran- 
domness for the 3D mode}2£i. In the inset the scaled distri- 
butions are shown as a function of y = (0 c (a, N) — 9 c )N llldv , 
see the text. 



In the second step we have studied the size-dependence 
of the distributions of the pseudo-critical points, which 
is illustrated in FigJ2]for the 3D model. From the scaling 
of the width, A6 C (N) ~ N~ 1 / d " w , and from the scaling 

tf-l/dv. we nave q]-)_ 

and v Sl respectively. 



of the mean value: \9 C — 8 C (N)\ ~ 
tained the critical exponents, v u 
We have calculated size-dependent effective exponents 
by two-point fits (comparing the results for sizes iV and 
N/2 d ), which are then extrapolated. The effective expo- 
nents for the 3D model are shown in the inset of FigO 
for the two different randomncsses. As in this example 
we have generally observed that the extrapolated critical 
exponents are universal, i.e. randomness independent. 
Estimates of the exponents are presented in Table [H to- 
gether with the values of the true critical points, 9 C 
and 9^J\ for the two randomnesses, respectively. One 
can notice in this Table that the error of the estimates 
is increasing with the dimensionality and the finite-size 
corrections are considerably strong for 4D and for ER 
random graphs. 

Having accurate estimates for the critical points we 
have renormalized the systems at 6 C and studied the 
scaling behavior of the moment of the correlation clus- 
ter, fj,(a, N), as well as that of the log-energy parameter, 
j(a,N) = —]ne(a,N). The average moment is found 
to scale as: Jl(N) ~ N d f' d , where df is the fractal di- 
mension of the correlation cluster. We illustrate this re- 
lation in Figj3] for the 3D and the 4D models, in which 
~p(N) is shown as a function of N in a log-log scale. In- 
deed, for not too small systems, N > 1000, the points 
are very well on straight lines, the slope of which be- 
ing the same for the two different randomnesses for the 
same d. From the cluster moment the magnetization is 
calculated as, m = ~p(N)/N, thus we have the scaling 
relations: x/d = j3/(dv) = 1 — df/d, where j3 is the 
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FIG. 3: (Color online) Average moment of the correlation 
cluster at the critical point vs. the size of the system in a log- 
log plot for the 3D and 4D models for two types of randomness 
(box-h: fixed-fe: +). The slope of the straight lines is 
given by: df/d = 0.387 and df/d = 0.32 for 3D and 4D, 
respectively. Inset: Finite-size effective exponents, v s (blue - 
dashed) and u w (red - full), for the 3D model for two types of 
randomness. In all cases the error of the calculation is smaller 
than the size of the symbol. 



magnetization exponent and x = (3 /v. Estimates for the 
exponents x/d, which are calculated through two-point 
fits are shown in Table HI 
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FIG. 4: (Color online) Distribution of the log-energy pa- 
rameters at the critical point in 3D, 4D and in the ER ran- 
dom graph for fixed-/i randomness and for different sizes (up- 
per panel)2£. In the lower panel the scaled distributions are 
shown, as described in the text. The constants in the scaling 
forms are 7o = -0.33 (3D), 70 = -0.23 (4D) and 70 = -0.5, 
lnTVo = -5. (ER). 

The energy parameter, e(a,N), is given by the value 
of the smallest effective transverse field, not considering 
that of the correlation cluster (if any). The distribution 
of the log-energy parameter, 7(0:, N), is shown in the up- 
per panel of FigUl for the different models. As a clear 
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TABLE I: Critical properties of the RTIM in three and four 
dimensions and in ER random graphs. N max denotes the 
number of spins in the largest finite systems used in the RG 
calculation. In the last four lines by estimating the differ- 
ent critical exponents results obtained by the two forms of 
disorder are taken into account. 





3D 


4D 


ER 


N max 


128 3 


48 4 


2 22 


op 


2.5305(10) 


3.110(5) 


2.775(2) 




-0.07627(2) 


-0.04698(10) 


-0.093(1) 


dv w 


2.90(15) 


3.30(15) 


7.(2) 


dv s 


2.96(5) 


2.96(15) 


5.(1) 


x/d 


0.613(5) 


0.68(3) 


0.81(4) 


iP 


0.46(2) 


0.46(2) 





indication of infinite disorder scaling the width of the 
distribution is increasing with N . In 3D and 4D the ap- 
propriate scaling variable is 7 = ('y(N) — j )N~^'/ d (j 
is a constant), as illustrated with the data collapse in 
the lower panel of FigfJJ The critical exponent ip nas 
been calculated from the optimal collapse of the dis- 
tributions, as well as from two-point fits by comparing 
the mean values 7~(iV) and 1 y(N/2 d ), which are presented 
in Table HI The ER random graphs arc infinite dimen- 
sional objects and in this case the broadening of the dis- 
tribution of the log-energy parameter is found to scale 
with In AT. As a good scaling combination we have here 
7 = (7(A0 - 7o)(ln iV/iVo)-" (70 and N are constants), 
which is illustrated with the data collapse in the lower 
panel of FigHJ with an exponent oj = 1.3(2). Thus the 
width of the distribution increases somewhat faster than 
linear in In AT, which fact justifies that also for ER ran- 
dom graphs the critical behavior of the RTIM is con- 
trolled by a logarithmically infinite disorder fixed point. 



IV. DISCUSSION 

Our numerical RG results indicate that the critical be- 
havior of the random transverse- field Ising model in three 
and four dimensions as well as in the ER random graph is 
controlled by infinite disorder fixed points. This fact jus- 
tifies the use of the SDRG method and ensures that the 
calculated numerical results about the critical exponents 
tend to be asymptotically correct for large sizes. Since 
the ER random graph represents the large-dimensional 
limit of the problem, infinite disorder scaling is expected 
to be valid at any dimensions. The critical exponents 
presented in Table Q] are found to be the same for the 
two types of ferromagnetic disorder used in our numeri- 
cal study. We expect therefore that the IDFP-s are at- 
tractive, at least for strong enough disorder. (For the 
weak-disorder behavior of the systems we can not make 
any definite statement using the SDRG method.) We 
note also that at the IDFP frustration does not matter, 
thus quantum spinglass and random quantum ferromag- 



net has the same infinite disorder fixed point. 

Singularities of the thermodynamic quantities at small 
temperatures involve the exponents in Table HI For ex- 
ample the susceptibility and the specific heat behave 
as: X (T) ~ (logT)^" 2 ^ /T and C V (T) ~ (logT) d /^, 
respectivel y 11 ' 24 . The analogous expressions for the ER 
random graph are: \og[T X (T)] ~ (1 - 2x/d)(log T) 1 ^ 
andlog[Cy(T)] ~ (logT) 1 /". 
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1/d 

FIG. 5: (Color online) Critical exponents of the RTIM as 
a function of 1/d. At 1/d = there are results of the ER 
random graph. 

As mentioned before the IDFP-s in Tabic Q] control the 
critical behavior of the spinglass transition, as well as 
that of a class of random quantum systems having an 
order parameter with discrete symmetry, such as ran- 
dom quantum Potts 3 -! and clock models^ 2 -. Noncquilib- 
rium phase transitions, such as the contact process with 
(strong) disorder— also belong to this class of universal- 
ity. The critical exponents in Table HI extending with the 
known results in 1D^ and 2Di&, show a smooth variation 
with the dimension, which is presented in FigJS] These 
results indicate that the large-o? limit of the problem is 
not singular. For a given d the correlation length criti- 
cal exponents, v s and u w , agree with each other, within 
the error of the calculation^ 3 .. These satisfy the rigor- 
ous bound 3 - 4 -, v > 2/d, and are in agreement with the 
scaling theory at conventional random fixed point o 35 ' 36 . 
Interestingly, the exponent tp is found very close to 1/2 
for any considered finite dimensional. This fact can be 
explained with our observation, that the low-energy ex- 
citations in any dimension are quasi-lD objects and the 
energy scale can be obtained by renormalizing these ob- 
jects practically independently of the rest of the system. 
This leads to approximately the same type of linear-size 
dependence of the energy in any dimension. 

At this point we want to mention a very recent study 
by Dimitrova and Mezard 3 ^ on the critical behavior of 
the RTIM by the cavity method^. In ID even the sim- 
ple mean-field cavity method is shown to recover some 
of the exact results, such as infinite disorder scaling at 
the critical point. On the contrary, results on the Bethe 
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lattice indicate the presence of a conventional random 
fixed point with a finite dynamical exponent. This result 
is probably due to the fact, that the local topology of 
the Bethe lattice is different from that of the hypercubic 
lattices, we considered in this paper. The local topology 
has already been found to have an important effect on 
the critical behavior of random quantum systems^. For 
example in the Bethe lattice it seems to be impossible 
to define an isotropic and quasi-one-dimcnsional cluster, 
which could be relevant for the low-energy excitations 
and thus for infinite disorder scaling. A direct SDRG 
study of the Bethe lattice RTIM could clarify some of 
the open questions. 

Finally we note, that the SDRG investigations pre- 
sented in this paper can be extended in several direc- 



tions. Here we mention the calculation of the entangle- 
ment entropy ^ 17 ' 41 in these systems, as well as study of 
the dynamical singularities in the disordered and ordered 
Griffiths phases^. 
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